#include <math.h>

float mutual_information(int length, int n1, int n2, int n12){
    /*
    Взаимная информация между двумя последовательностями длины l
    n1 - число "1" в 1й последовательности
    n2 - число "1" в 2й последовательности
    n12 - число "1" в пересечении последовательностей
    */

    float I;
    float l = float(length);
    float n12_th = n1 * n2 / float(l);
    float n0_th  = float(n1 - n12_th + n2 - n12_th);
    float n1_th  = l - n0_th;
    float p0_th = n0_th / l;
    float p1_th = n1_th / l;
    float n0_exp = float(n1 - n12 + n2 - n12);
    float n1_exp = l - n0_exp;
    if ((n0_exp > 0) and (n1_exp > 0) and (n0_th) > 0 and (n1_th > 0)) {
	I = n1_exp * log(n1_exp) + n0_exp * log(n0_exp) \
	    - n1_exp * log(p1_th) - n0_exp * log(p0_th) - l * log(l);
    } else {
	I = 0.0;
    }
    return I;
}

